1 Чтение и предобработка данных

Подгружаем датасет

led <- read_rds("./life_expectancy_data.RDS")

Смотрим на переменные

str(led)
## Classes 'data.table' and 'data.frame':   195 obs. of  23 variables:
##  $ Country                                : chr  "Afghanistan" "Albania" "Algeria" "Angola" ...
##  $ Year                                   : int  2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
##  $ Gender                                 : chr  "Female" "Female" "Female" "Female" ...
##  $ Life expectancy                        : num  66.4 80.2 78.1 64 78.1 ...
##  $ Unemployment                           : num  14.06 11.32 18.63 7.84 8.26 ...
##  $ Infant Mortality                       : num  42.9 7.7 18.6 44.5 5.1 ...
##  $ GDP                                    : num  18799450743 15400242875 171767403748 89417190341 1687533333 ...
##  $ GNI                                    : num  19098305631 15198661275 167568133680 81900890781 1581477852 ...
##  $ Clean fuels and cooking technologies   : num  36 80.7 99.3 49.6 100 ...
##  $ Per Capita                             : num  494 5396 3990 2810 17377 ...
##  $ Mortality caused by road traffic injury: num  15.9 11.7 20.9 26.1 0 ...
##  $ Tuberculosis Incidence                 : num  189 16 61 351 0 29 26 2.2 6.9 6 ...
##  $ DPT Immunization                       : num  66 99 91 57 95 ...
##  $ HepB3 Immunization                     : num  66 99 91 53 99 ...
##  $ Measles Immunization                   : num  64 95 80 51 93 ...
##  $ Hospital beds                          : num  0.432 3.052 1.8 0.8 2.581 ...
##  $ Basic sanitation services              : num  49 99.2 86.1 51.4 85.5 ...
##  $ Tuberculosis treatment                 : num  91 88 86 69 72.3 ...
##  $ Urban population                       : num  25.8 61.2 73.2 66.2 24.5 ...
##  $ Rural population                       : num  74.2 38.8 26.8 33.8 75.5 ...
##  $ Non-communicable Mortality             : num  36.2 6 12.8 19.4 17.6 ...
##  $ Sucide Rate                            : num  3.6 2.7 1.8 2.3 0.8 ...
##  $ continent                              : Factor w/ 5 levels "Africa","Americas",..: 3 4 1 1 2 2 4 2 5 4 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "Country"
summary(led)
##    Country               Year         Gender          Life expectancy
##  Length:195         Min.   :2019   Length:195         Min.   :55.49  
##  Class :character   1st Qu.:2019   Class :character   1st Qu.:70.02  
##  Mode  :character   Median :2019   Mode  :character   Median :77.55  
##                     Mean   :2019                      Mean   :75.52  
##                     3rd Qu.:2019                      3rd Qu.:80.95  
##                     Max.   :2019                      Max.   :88.10  
##   Unemployment    Infant Mortality      GDP                
##  Min.   : 0.178   Min.   : 1.40    Min.   :     188391771  
##  1st Qu.: 3.735   1st Qu.: 5.35    1st Qu.:   11170154590  
##  Median : 5.960   Median :13.50    Median :   39670977333  
##  Mean   : 8.597   Mean   :19.61    Mean   :  465995962784  
##  3rd Qu.:10.958   3rd Qu.:30.23    3rd Qu.:  247555140295  
##  Max.   :36.442   Max.   :75.80    Max.   :21433224697000  
##       GNI                 Clean fuels and cooking technologies
##  Min.   :     375392118   Min.   :  0.00                      
##  1st Qu.:   10943629972   1st Qu.: 34.50                      
##  Median :   40087702395   Median : 80.70                      
##  Mean   :  486438131433   Mean   : 65.98                      
##  3rd Qu.:  245717112337   3rd Qu.:100.00                      
##  Max.   :21708650000000   Max.   :100.00                      
##    Per Capita       Mortality caused by road traffic injury
##  Min.   :   228.2   Min.   : 0.00                          
##  1st Qu.:  2165.3   1st Qu.: 8.20                          
##  Median :  6624.8   Median :16.00                          
##  Mean   : 16821.0   Mean   :17.06                          
##  3rd Qu.: 19439.7   3rd Qu.:24.00                          
##  Max.   :175813.9   Max.   :64.60                          
##  Tuberculosis Incidence DPT Immunization HepB3 Immunization
##  Min.   :  0.0          Min.   :35.00    Min.   :35.00     
##  1st Qu.: 12.0          1st Qu.:85.69    1st Qu.:81.31     
##  Median : 46.0          Median :92.00    Median :91.00     
##  Mean   :103.8          Mean   :87.99    Mean   :86.76     
##  3rd Qu.:138.5          3rd Qu.:97.00    3rd Qu.:96.00     
##  Max.   :654.0          Max.   :99.00    Max.   :99.00     
##  Measles Immunization Hospital beds    Basic sanitation services
##  Min.   :37.00        Min.   : 0.200   Min.   :  8.632          
##  1st Qu.:84.85        1st Qu.: 1.301   1st Qu.: 62.919          
##  Median :92.00        Median : 2.570   Median : 91.144          
##  Mean   :87.31        Mean   : 2.997   Mean   : 77.380          
##  3rd Qu.:96.50        3rd Qu.: 3.773   3rd Qu.: 98.582          
##  Max.   :99.00        Max.   :13.710   Max.   :100.000          
##  Tuberculosis treatment Urban population Rural population
##  Min.   :  0.00         Min.   : 13.25   Min.   : 0.00   
##  1st Qu.: 73.00         1st Qu.: 41.92   1st Qu.:21.98   
##  Median : 82.00         Median : 58.76   Median :41.24   
##  Mean   : 77.57         Mean   : 59.12   Mean   :40.88   
##  3rd Qu.: 88.00         3rd Qu.: 78.02   3rd Qu.:58.08   
##  Max.   :100.00         Max.   :100.00   Max.   :86.75   
##  Non-communicable Mortality  Sucide Rate        continent 
##  Min.   : 4.40              Min.   : 0.300   Africa  :52  
##  1st Qu.:11.85              1st Qu.: 2.050   Americas:38  
##  Median :17.20              Median : 3.500   Asia    :42  
##  Mean   :17.05              Mean   : 4.802   Europe  :48  
##  3rd Qu.:22.10              3rd Qu.: 6.600   Oceania :15  
##  Max.   :43.70              Max.   :30.100

Все выглядит пристойно.

Оценим пропущенные значения.

led %>% 
  is.na() %>% 
  sum()
## [1] 0

Их нет, и это прекрасно.

EDA

2

Посмотрим на взаимное распределение встречаемости туберкулеза и частоты самоубийств. Добавил цветом принадлежность государства к одному из континентов.

# Задаем палетку:
pal <- c("#DF536B", "#61D04F", "#2297E6", "#28E2E5", "#CD0BBC")
pal <- setNames(pal, levels(led$continent))

# Рисуем график
plot_ly(
  led,
  x = ~ `Tuberculosis Incidence`,
  y = ~ `Sucide Rate`,
  size = 1,
  color = ~ continent,
  colors = pal,
  alpha = 0.5,
  width = 1.5,
  height = 1.5
) %>% 
  layout(
      title = "Распределение частоты самоубийств в зависимости от\nвстречаемости туберкулеза у жителей различных континентов",
      xaxis = list(title = "Встречаемость туберкулеза",
                   zeroline = FALSE),
      yaxis = list(title = "Частота самоубийств",
                   zeroline = FALSE)
  )
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Может, рисунок не очень хорош с точки зрения анализа данных, но как картина в духе пуантилизма он вполне неплох.

Анализ данных

3

Проверим, различаются ли средняя ожидаемая продолжительность жизни в странах Африки и Америки.

Для сравнения применим t-тест Стьюдента. Чтобы понять, необходимо ли нам применять тест Уэлча, оценим равенство дисперсий.

leveneTest(`Life expectancy` ~ continent, data = led %>% filter(continent %in% c("Africa", "Americas")))
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  1  7.5392 0.007319 **
##       88                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дисперсии негомогенны. Применим t-тест с поправкой Уэлча (t-тест по умолчанию в R)

for_vis_tt <- led %>% 
  filter(continent %in% c("Africa", "Americas")) %>%
  rstatix::t_test(`Life expectancy` ~ continent) %>% 
  add_xy_position(x = "continent")

И визуализируем его:

led %>% filter(continent %in% c("Africa", "Americas")) %>% 
ggplot(aes(x=continent, y = `Life expectancy`))+
  geom_boxplot(fill="green")+
  labs(
    x = "Континент",
    y = "Ожидаемая\nпродолжительность жизни",
    title = "Распределение ожидаемой продолжительности жизни\nв странах Африки и Америки"
  )+
  theme_bw()+
  theme(
    axis.text.x = element_text(size=16),
    axis.text.y = element_text(size=16),
    axis.title.x = element_text(size=19), 
    axis.title.y = element_text(size=19),
    plot.title = element_text(size=23, hjust=5),
    legend.title = element_text(size=19),
    legend.text = element_text(size=16),
    legend.position = "bottom",
    legend.key.size = unit(1.5, "line"),
    strip.text.x = element_text(size=16)
  )+
  stat_pvalue_manual(for_vis_tt, tip.length = 0)

Средняя ожидаемая продолжительность жизни в странах Африки и Америки различаются.

4

Сделаем корреляционный анализ числовых переменных.

led_num <- led %>% 
  select(where(is.numeric) & !Year)

ggpairs(led_num, progress = FALSE)

График для 19-ти переменных выглядит довольно громоздко. Построим коррплот.

# Спирменовская корреляция -- потому что мы не проверяли распределение данных, а наверняка там найдется что-то ненормально-нелинейное
led_num_corr <- corr.test(led_num , method="spearman")
corr_plot <- ggcorrplot(led_num_corr$r, 
                        type = "lower", 
                        outline.col = "white", 
                        lab=TRUE, 
                        method = "circle", 
                        p.mat = led_num_corr$p,
                        title = "Скоррелированность количественных данных\nв датасете по ожидаемой продолжительности жизни",
                        legend.title = "Коэффициент\nкорреляции")
corr_plot

Наблюдается большое количество корреляций. Для примера можно привести сильную положительную скоррелированность между собой показателей вакцинации; продолжительности жизни с подушевым доходом и санитарно-гигиеническими нормами. Самая сильная (обратная) корреляция, что логично, наблюдается между урбанизацией и долей сельского населения.

5

Ростроим иерархическую кластеризацию

# Делаем матрицу дистанций
led_num_scaled <- scale(led_num)

led_num_dist <- dist(led_num_scaled, 
                        method = "euclidean")
led_num_hc <- hclust(d = led_num_dist, 
                        method = "ward.D2")
fviz_dend(led_num_hc, 
          cex = 0.1)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

6

Реализуем фитмэп

pheatmap(led_num_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = led_num_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5,   # Прикидочное значение
         cutree_cols = length(colnames(led_num_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

Мы наблюдаем следующие 2 группы скоррелированных переменных: 1) наличие иммунизации от кори, гепатита B и АКДС 2) экономические метрики: ВВП и ВНД

И 2 группы менее скоррелированных: 1) Безработица, встречаемость туберкулеза, лечение туберкулеза, смертность в ДТП, детская смертность и доля сельского населения. 2) Частота самоубийств, количество больничных коек (страшная вырисовывается корреляция), базовые санитарарные услуги, продолжительность жизни, доля городского населения, подушевой доход.

В данных мы выделили 4 субкластера (некий гиперпараметр), один из которых чрезвычайно мал (второй сверху на вертикальной дендрограмме) и характеризуется высокими значениями экономических показателей. Первый субкластер отличается очень низкими значениями иммунизации. Остальные три кластера выглядят достаточно негомогенными по большинству признаков (кроме экономических показателей).

7

Построим PCA

led_full.pca <- prcomp(led_num_scaled) # Это уже отшкалированные данные, дополнительно не шкалируем

summary(led_full.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18                  PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 0.0000000000000002735
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.0000000000000000000
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.0000000000000000000

Первые 5 компонент объясняют 75% вариации. Не идеально, но приемлемо. Посмотрим на соответствующий график.

fviz_eig(led_full.pca, addlabels = T, ylim = c(0, 40))

Посмотрим на переменные

fviz_pca_var(led_full.pca, col.var = "contrib")

Заметим, что все так же выделяется группа иммунизации, и две группы “слабо скоррелированных” переменных с фитмэпа.

fviz_pca_var(led_full.pca, 
             select.var = list(contrib = 5),
             col.var = "contrib")

5 переменных, вносящих наибольший вклад – все те же три переменные о проценте иммунизации, продолжительность жизни и младенческая смертность.

Посмотрим, кто вкладывается в избранные 6 главных компонент

fviz_contrib(led_full.pca, choice = "var", axes = 1, top = 24) # 1

fviz_contrib(led_full.pca, choice = "var", axes = 2, top = 24) # 2

fviz_contrib(led_full.pca, choice = "var", axes = 3, top = 24) # 3

fviz_contrib(led_full.pca, choice = "var", axes = 4, top = 24) # 4

fviz_contrib(led_full.pca, choice = "var", axes = 5, top = 24) # 5

В первую наибольший вклад вносят продолжительность жизни, младенческая смертность, санусловия и соблюдение гигиены и правил приготовления пищи.

Во вторую – переменные вакцинации.

В третью – подавляющий вклад со стороны экономических показателей.

В четвертую – частота самоубийств.

В пятую основной вклад вносит безработица.

8

Построим би-2плот

pca_plot <- ggbiplot::ggbiplot(led_full.pca, 
         scale=0, 
         groups = as.factor(led$continent), 
         labels = led$Country,
         labels.size = 2,
         alpha = 0.5,
         varname.size = 3,
         varname.abbrev = TRUE,
         varname.adjust = 0) + 
  theme_minimal()
ggplotly(pca_plot,
         tooltip = "label")

9

В ходе анализа биплота выявляются следующие закономерности: 1. В европейских странах иммунизация проводится с большой частотой. В африканских – увы. 2. Азиатские страны представляют собой гетерогенную во всех отношениях группу. 3. Субкластер американских стран, Китай и его автономии имеют высокий уровень урбанизации и экономиечксих показателей (США и вовсе выглядят выбросом). 4. Большая часть африканских стран имеет высокую встречаемость туберкулеза, смертности от неинфекционных заболеваний, долю сельского населения и младенческую смертность. Это сопровождается более низкой ожидаемой продолжительностью жизни, меньшим числом больничных коек, более низким показателем чистоты топлива и технологий для приготовления пищи. 5. Наибольшая ожидаемая продолжительность жизни – в Японии. 6. Если вы не ищете в жизни легких путей – добро пожаловать в Сомали.

10

Построим UMAP на тех же данных.

led_umap <- led %>% 
  select(where(is.numeric) | continent | Country)

umap_prep <- recipe(~., data = led_num) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors()) %>%  
  prep() %>%
  juice()
um <- umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(colour = led_umap$continent), 
             alpha = 0.7, size = 2) +
  labs(color = NULL) 
ggplotly(um)

На UMAP мы наблюдаем два облака точек. Одно включает в себя большинство африканских стран и несколько стран остальных регионов (кроме Европы).

Второе включает в себя все европейские страны, большую часть американских и азиатских, а также несколько стран Океании и Африки.

Отличия графиков.

  1. UMAP в отличие от PCA без фиксации seed невоспроизводим. Каждый раз он принимает несколько иную конформацию (иногда появляется третий маленький промежуточный и гетерогенный кластер).
  2. Визуализация UMAP сильно зависит от гиперпарметров (числа соседей и прочего), PCA же константен на одних и тех же данных.
  3. Это обуславливает тот аспект, что PCA может применять для анализа и построения умозаключений, а UMAP перимущественно для EDA.
  4. Визуально на UMAP кластера достаточно отдалены друг от друга (если, конечно, не задать им специфические настройки) за счет того, что этот метод в первую очередь ориентируется на расстояние между соседними точками, а не на глобальную картину в целом. PCA выглядит более диффузно. Говоря образно, можно сравнить UMAP со ступенчатым переходом, а PCA – с континуумом.

11

Будем случайно удалять 5 переменных и строить PCA.

Попытка 1

Трансформируем данные

set.seed(1929)
led_rand_1 <- led_num %>% 
  select(sample(1:ncol(led_num), size=ncol(led_num)-5))
led_rand_1_scaled <- scale(led_rand_1)
led_rand_1.pca <- prcomp(led_rand_1_scaled) # Это уже отшкалированные данные, дополнительно не шкалируем

summary(led_rand_1.pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.5004 1.4544 1.12167 1.00271 0.91350 0.82904 0.76124
## Proportion of Variance 0.4466 0.1511 0.08987 0.07182 0.05961 0.04909 0.04139
## Cumulative Proportion  0.4466 0.5977 0.68753 0.75935 0.81895 0.86805 0.90944
##                            PC8     PC9   PC10   PC11    PC12    PC13
## Standard deviation     0.68694 0.59450 0.4267 0.3588 0.29900 0.20587
## Proportion of Variance 0.03371 0.02525 0.0130 0.0092 0.00639 0.00303
## Cumulative Proportion  0.94314 0.96839 0.9814 0.9906 0.99697 1.00000
##                                         PC14
## Standard deviation     0.0000000000000002611
## Proportion of Variance 0.0000000000000000000
## Cumulative Proportion  1.0000000000000000000
fviz_eig(led_rand_1.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(led_rand_1.pca, col.var = "contrib")

pca_plot <- ggbiplot::ggbiplot(led_rand_1.pca, 
         scale=0, 
         groups = as.factor(led$continent), 
         labels = led$Country,
         labels.size = 2,
         alpha = 0.5,
         varname.size = 3,
         varname.abbrev = TRUE,
         varname.adjust = 0) + 
  theme_minimal()
ggplotly(pca_plot,
         tooltip = "label")

Попытка 2

set.seed(1949)
led_rand_2 <- led_num %>% 
  select(sample(1:ncol(led_num), size=ncol(led_num)-5))
led_rand_2_scaled <- scale(led_rand_2)
led_rand_2.pca <- prcomp(led_rand_2_scaled) # Это уже отшкалированные данные, дополнительно не шкалируем

summary(led_rand_2.pca)
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.532 1.3391 1.16256 1.03700 0.95258 0.86420 0.69783
## Proportion of Variance 0.458 0.1281 0.09654 0.07681 0.06482 0.05335 0.03478
## Cumulative Proportion  0.458 0.5860 0.68259 0.75940 0.82422 0.87756 0.91235
##                            PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.62484 0.59731 0.46624 0.38057 0.27657 0.20309
## Proportion of Variance 0.02789 0.02548 0.01553 0.01035 0.00546 0.00295
## Cumulative Proportion  0.94023 0.96572 0.98124 0.99159 0.99705 1.00000
##                                         PC14
## Standard deviation     0.0000000000000002401
## Proportion of Variance 0.0000000000000000000
## Cumulative Proportion  1.0000000000000000000
fviz_eig(led_rand_2.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(led_rand_2.pca, col.var = "contrib")

pca_plot <- ggbiplot::ggbiplot(led_rand_2.pca, 
         scale=0, 
         groups = as.factor(led$continent), 
         labels = led$Country,
         labels.size = 2,
         alpha = 0.5,
         varname.size = 3,
         varname.abbrev = TRUE,
         varname.adjust = 0) + 
  theme_minimal()
ggplotly(pca_plot,
         tooltip = "label")

Попытка 3

set.seed(1988)
led_rand_3 <- led_num %>% 
  select(sample(1:ncol(led_num), size=ncol(led_num)-5))
led_rand_3_scaled <- scale(led_rand_3)
led_rand_3.pca <- prcomp(led_rand_3_scaled) # Это уже отшкалированные данные, дополнительно не шкалируем

summary(led_rand_3.pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3912 1.4670 1.15381 1.05385 0.97230 0.90148 0.81113
## Proportion of Variance 0.4084 0.1537 0.09509 0.07933 0.06753 0.05805 0.04699
## Cumulative Proportion  0.4084 0.5622 0.65724 0.73657 0.80410 0.86215 0.90914
##                            PC8     PC9    PC10    PC11    PC12    PC13
## Standard deviation     0.67221 0.57069 0.47125 0.37708 0.29838 0.20289
## Proportion of Variance 0.03228 0.02326 0.01586 0.01016 0.00636 0.00294
## Cumulative Proportion  0.94142 0.96468 0.98054 0.99070 0.99706 1.00000
##                                         PC14
## Standard deviation     0.0000000000000003701
## Proportion of Variance 0.0000000000000000000
## Cumulative Proportion  1.0000000000000000000
fviz_eig(led_rand_3.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(led_rand_3.pca, col.var = "contrib")

pca_plot <- ggbiplot::ggbiplot(led_full.pca, 
         scale=0, 
         groups = as.factor(led$continent), 
         labels = led$Country,
         labels.size = 2,
         alpha = 0.5,
         varname.size = 3,
         varname.abbrev = TRUE,
         varname.adjust = 0) + 
  theme_minimal()
ggplotly(pca_plot,
         tooltip = "label")

В первых двух случаях мы достигаем порога объясненной дисперсии в 75% на 4-й компоненте. В 3-м – между 4-й (73%) и 5-й (80%). “Стрелка-плоты” во всех трех случаях отличаются. Во втором случае сложнее сформировать ярко выраженные группы переменных (оттого и точки распределены более дисперсно). В первом же мы имеем иное направление компонент, из-за чего график несколько меняет свою “кранио-каудальную” структуру. Вместе с тем, основные паттерны полной версии PCA на них сохраняются (но, например, исчезает США-выброс). Третий вариант наиболее близок к оригинальной версии PCA.

Исключение отдельных переменных меняет вклад в общую дисперсию и структуру скоррелированности данных. Выпадение важных для объяснения структуры данных переменных может резко переменить отображение PCA.

12

Добавляем колокни типа: да-нет

led_dummy <- led %>% 
  select(where(is.numeric) | continent & !Year) %>% 
  mutate(oceania = if_else(continent == "Oceania", 1, 0),
         africa = if_else(continent == "Africa", 1, 0)
  ) %>% 
  select(-c(continent, Year)) %>% 
  scale
led_dummy.pca <- prcomp(led_dummy) # Это уже отшкалированные данные, дополнительно не шкалируем

summary(led_dummy.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion  0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion  0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
##                           PC15    PC16    PC17    PC18    PC19    PC20
## Standard deviation     0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion  0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
##                                         PC21
## Standard deviation     0.0000000000000002623
## Proportion of Variance 0.0000000000000000000
## Cumulative Proportion  1.0000000000000000000
fviz_eig(led_dummy.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(led_dummy.pca, col.var = "contrib")

pca_plot <- ggbiplot::ggbiplot(led_full.pca, 
         scale=0, 
         groups = as.factor(led$continent), 
         labels = led$Country,
         labels.size = 2,
         alpha = 0.5,
         varname.size = 3,
         varname.abbrev = TRUE,
         varname.adjust = 0) + 
  theme_minimal()
ggplotly(pca_plot,
         tooltip = "label")

Введение дамми-переменных выглядит нерациональным при текущем подходе. Они представлены в бинарном числовом формате, а после шкалирования превращаются в нечто совершенно невообразимое. Они не объясняют дополнительно часть дисперсии данных и не очень хорошо скоррелированы с остальными параметрами (что видно на биплоте и стрелка-плоте). Они не покрывают все 5 возможных вариантов континентов.